home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / poly / solve_cubic.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  2.9 KB  |  111 lines
  1. /* poly/solve_cubic.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* solve_cubic.c - finds the real roots of x^3 + a x^2 + b x + c = 0 */
  21.  
  22. #include <config.h>
  23. #include <math.h>
  24. #include <gsl/gsl_math.h>
  25. #include <gsl/gsl_poly.h>
  26.  
  27. #define SWAP(a,b) do { double tmp = b ; b = a ; a = tmp ; } while(0)
  28.  
  29. int 
  30. gsl_poly_solve_cubic (double a, double b, double c, 
  31.                       double *x0, double *x1, double *x2)
  32. {
  33.   double q = (a * a - 3 * b);
  34.   double r = (2 * a * a * a - 9 * a * b + 27 * c);
  35.  
  36.   double Q = q / 9;
  37.   double R = r / 54;
  38.  
  39.   double Q3 = Q * Q * Q;
  40.   double R2 = R * R;
  41.  
  42.   double CR2 = 729 * r * r;
  43.   double CQ3 = 2916 * q * q * q;
  44.  
  45.   if (R == 0 && Q == 0)
  46.     {
  47.       *x0 = - a / 3 ;
  48.       *x1 = - a / 3 ;
  49.       *x2 = - a / 3 ;
  50.       return 3 ;
  51.     }
  52.   else if (CR2 == CQ3) 
  53.     {
  54.       /* this test is actually R2 == Q3, written in a form suitable
  55.          for exact computation with integers */
  56.  
  57.       /* Due to finite precision some double roots may be missed, and
  58.      considered to be a pair of complex roots z = x +/- epsilon i
  59.      close to the real axis. */
  60.  
  61.       double sqrtQ = sqrt (Q);
  62.  
  63.       if (R > 0)
  64.     {
  65.       *x0 = -2 * sqrtQ  - a / 3;
  66.       *x1 = sqrtQ - a / 3;
  67.       *x2 = sqrtQ - a / 3;
  68.     }
  69.       else
  70.     {
  71.       *x0 = - sqrtQ  - a / 3;
  72.       *x1 = - sqrtQ - a / 3;
  73.       *x2 = 2 * sqrtQ - a / 3;
  74.     }
  75.       return 3 ;
  76.     }
  77.   else if (CR2 < CQ3) /* equivalent to R2 < Q3 */
  78.     {
  79.       double sqrtQ = sqrt (Q);
  80.       double sqrtQ3 = sqrtQ * sqrtQ * sqrtQ;
  81.       double theta = acos (R / sqrtQ3);
  82.       double norm = -2 * sqrtQ;
  83.       *x0 = norm * cos (theta / 3) - a / 3;
  84.       *x1 = norm * cos ((theta + 2.0 * M_PI) / 3) - a / 3;
  85.       *x2 = norm * cos ((theta - 2.0 * M_PI) / 3) - a / 3;
  86.       
  87.       /* Sort *x0, *x1, *x2 into increasing order */
  88.  
  89.       if (*x0 > *x1)
  90.     SWAP(*x0, *x1) ;
  91.       
  92.       if (*x1 > *x2)
  93.     {
  94.       SWAP(*x1, *x2) ;
  95.       
  96.       if (*x0 > *x1)
  97.         SWAP(*x0, *x1) ;
  98.     }
  99.       
  100.       return 3;
  101.     }
  102.   else
  103.     {
  104.       double sgnR = (R >= 0 ? 1 : -1);
  105.       double A = -sgnR * pow (fabs (R) + sqrt (R2 - Q3), 1.0/3.0);
  106.       double B = Q / A ;
  107.       *x0 = A + B - a / 3;
  108.       return 1;
  109.     }
  110. }
  111.